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The above pictures show the analytical and computational aspects of turbulence and combustion. These figures show an un- 
stable detonation wave in two distinct states. The top figure shows the temperature profile of a Z-N-D detonation wave with qo 
(the heat release parameter) of 50, E* — activation energy) of 10, and f (the overdriver of the detonation) of 1.2. Note the 


eneration of regularly paired eddies. The bottom picture is produced under the same parameter values except for E* that 
as a value of 50. This picture shows the generation of strong irregular cells and fully developed 2-D turbulence. 


This pioneering research of Professor Andrew Majda and co-workers at Princeton University was sponsored by the Office of 
Naval Research. On page 20, Professor Majda is featured in “Profiles in Science.” 


ONR has been working in this field of research for 30 years. The discovery in the 1960's under ONR support of coherence in 
turbulence revolutionized the view of turbulent flow. Before that time the statistical description of turbulence emphasized its 
randomness and brought about time-averaged turbulence models. ONR-sponsored research identified coherent structures in 
both free and bounded shear flows, and subsequent programs investigated the role of these structures in the production and 
dissipation of energy in the flow. 
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Increasing, mixing, and generating the optimum large scale structures and small scale turbulence are crucial to increase 
the combustion efficiency of ramjets (photograph on the left hand side of the cover). Non axisymmetric nozzles have 
been found to generate enhanced mixing by the phenomena called “axis switching,” where the vortex rings deform to 
shift axes (major axis to minor axis). This causes the jet to spread more, generating increased mixing. On the right hand 
side is shown a square jet with the deformed vortex rings at the bottom, and the longitudinal vortices and their merging 
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Research on Multiple- 
Valued Logic At The 
Naval Postgraduate 


School 


Jon T. Butler 


Department of Electrical and Computer Engineering 


Naval Postgraduate School 


Abstract 


The Navy’s need for faster, smaller computers has in- 
spired research on the use of more than two levels of logic in 
computer circuits. This paper focusses on work at NPS in 
multiple-valued logic. It discusses the development of multi- 
ple-valued programmable logic arrays (MVL-PLA), as well 
as computer-aided design tools for MVL-PLA’s. 


Introduction 


The effectiveness of Navy systems is dependent on com- 
puters. Computers are used to predict weather, to track re- 
sources, to guide missiles, to simulate combat, and to process 
the words and numbers necessary for the general conduct of 
Navy business. Their usefulness depends on computing 
power, speed, and size. As with other critical resources, the 
Navy cannot depend on the private sector as the exclusive 
provider of the technology; it must be an active participant. 
This article focuses on computer circuits, and specifically on 
multiple-valued circuits, whose use improves computing 
power, speed, and size. We first consider binary circuits. Then, 
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we discuss multiple-valued circuits. This serves as introduc- 
tion to the main discussion, research at the Naval Postgraduate 
School on this important subject. 


Disadvantages of Binary 


Current computers are binary; they use two values of 
logic. For convenience, we represent these as 0 and 1. But 0 
and 1 are simply symbols for some quantity, such as voltage, 
current, or charge. For example, in the most widely used 
circuits, 0 represents 0 volts and 1 represents 5 volts. 

Binary is commonly used for historical reasons. Comput- 
ers in the 1940’s used relays, which had two stable states, open 
and closed. Tubes and transistors have two stable states, satu- 
ration (conducting) and cutoff (nonconducting). There are dis- 
advantages to binary, however. One is evident by the number of 
bits needed to represent a number. For example, the decimal 
number 33 corresponds to 10001 in binary. That is, while five 
symbols are needed to represent 33 in binary, only two are needed 
in decimal. Unfortunately, the difference increases with the size 
of the number. This inherent complexity is hidden from most 
users by hardware (e.g., binary to decimal converters in digital 





watches) and by software (e.g., high level languages). How- 
ever, certain disadvantages cannot be so simply hidden. Real 
penalties exist in both compactness and speed. 


Compactness 


The complexity of the circuits that perform addition, 
multiplication, and other arithmetic operations depends on the 
number of symbols used. For commonly used number sys- 
tems, this complexity increases with the number of symbols. 
As aresult, so does the chip area needed for devices. However, 
the most significant penalty is not with the devices, but the 
interconnect between devices. In VLSI circuits, about 70% of 
the chip area is used for interconnect, 20% for insulation, and 
10% for devices. Interconnect corresponds to the chip area 
used to carry information around the circuit. Insulation is 
needed to separate devices and interconnect. With so much 
area used to connect devices, less area is available for logic. 

The inefficient use of interconnect by binary also occurs 
outside the chip. Signals going into and out of a VLSI chip 
must use pins connecting the chip to the circuit board. This is 
called the pinout problem. That is, while significant reductions 
have occurred in the size of logic devices, there has been no 
comparable reduction in the size of integrated circuit pin 
connections. Strength and reliability dictate a (relatively large) 
minimum pin connection size. As a result, there is a need to 
better use existing pins. 


Speed 


A limit on the speed of arithmetic circuits is the carry (or 
borrow) between digits. For example, in binary addition, the 
most significant sum bit depends not only on the most signif- 
icant bits of the two numbers being added, but also on all lower 
bits (contrast this to the least significant sum bit which does 
not depend on higher order bits). The dependence occurs 
through the carry between bits. Thus, computation of the most 
significant sum bit must wait until the carries out of all lower 
order bits are computed. Addition, when done in this way, 
cannot be performed at high speed. 

However, there are multiple-valued number systems 
where this dependence does not exist. An example is the 
residue number system. In this system, operations that form 
the sum, difference, or product occur only at each digit inde- 
pendent of all other digits. Because of this independence, 
arithmetic operations are fast in the 1esidue number system. 


Multiple-Valued Logic 


All of the disadvantages discussed above can be alleviated 
by using more than two levels of logic. The problem of a large 
number of bits needed to represent a number decreases as the 
number of logic values increases. Interconnect, both internal and 
external to the VLSI chip, is more efficiently used in multiple- 





Figure 1. 


Photomicrograph of a Basic Charge Coupled Device Program- 
mable Logic Array (Courtesy of the University of Twente, OIEEE 
1988). 








valued logic. The carry propagation problem disappears with 
the proper choice of multiple-valued number system. 

However, there remains the question of implementing a 
multiple-valued system. Modern computing devices naturally 
allow more than two levels of logic. For example, in charge- 
coupled devices (CCD), logic levels are represented as various 
quantities of charge”. Resonant-tunnelling devices? allow four 
or more levels of voltage. Fuzzy logic devices offer large 
numbers of logic levels; so many, that they are modelled as 
infinite; present circuit implementations exceed 15 levels’’. 
Bio-molecular devices, where information processing occurs 
at the molecule level®, operate on an enormous number of logic 
levels, conservatively estimated to be over 10,000. 


Multiple-Valued Programmable Logic Arrays 


The design of logical circuits is a complex process. The 
vast amount of circuitry needed for even a simple operation 
like addition can be daunting. Because of this, a way is needed 
to organize circuits. The programmable logic array or PLA is 
a Circuit structure designed to handle this complexity. PLA’s 
are uniform and thus easily replicated in VLSI. They can be 
programmed, and therefore adapted to realize a given design 
specification. Because of this characteristic, our research at the 
Naval Postgraduate School has concentrated on PLA design. 

In a cooperative venture with the University of Twente in 
Enschede, The Netherlands, we have developed a PLA circuit 
for charge-coupled device (CCD) technology (See Figure 1 
from’). In CCD, values are stored in cells (capacitors) as 
charge. A logic 0, 1, 2, and 3 is stored as 0, 1,000,000, 
2,000,000, and 3,000,000 electrons. Computation occurs as 
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charge moves among cells combining and breaking up. An 
advantage of CCD is compactness, more so than any other 
VLSI technology. While slower than the common CMOS 
technology, it is faster than disk and stands as a replacement 
for disk. The use of multiple-valued logic in CCD increases 
its density significantly. CCD is used as a light receptor, and 
the advent of logic circuits in CCD means that image process- 
ing can be done on the same chip in which the image is 
received. With large numbers of logic levels, there is no need 
for encoders between light receptors and the logic used to 
process and store the image. CCD circuits with many logic 
levels have been fabricated. For example, Hitachi has suc- 
ceeded in building a 16 level CCD memory. 

While the uniform structure of a PLA helps one to handle 
the high complexity of logic circuits, there is the problem of 
designing large PLA’s. To make multiple-valued PLA design 
feasible for large systems, we initiated a program to develop 
computer-aided design tools. 


Computer-Aided Design for 
Multiple-Valued Logic 


We began this program with a study of heuristic tech- 
niques for finding minimal PLA’s, PLA’s with the smallest 
structure. A heuristic will produce a design, although not 
necessarily a minimal one. In 1987 at the beginning of this 
Study, there were three heuristic techniques. At that time, there 
was no systematic analysis of their relative merits; each ex- 
isted separately. In an effort to gain an understanding of the 
multiple-valued PLA design process, we compared the three 
heuristics on the same functions (specifications). One interest- 
ing result’? was that heuristics tended to perform better on 
different classes of functions, and there was a surprising 
advantage to applying all three and choosing the most compact 
realization. 

We also provided a partial answer to an important ques- 
tion; How well do heuristics perform compared to the best that 
one could possibly do? To answer this, we developed the first 
practical approach to finding provably minimal solutions; an 
efficient search technique that is much faster than exhaustive 
search, but still guaranteeing a minimal solution. Indeed, 
without the improved search technique, we would have been 
unable to achieve the minimal solution because of the enor- 
mous computation time required by exhaustive search. This 
comparison provided much needed insight into heuristic meth- 
ods and provided realism in our expectation of the quality of 
solutions provided by heuristic methods. 

It is important to state that these experiments deal with 
small functions; for even medium size functions this experi- 
ment would have been impossible. For larger functions, the 
only approach is heuristic. Indeed, our research shows that 
often a minimal realization is not produced. Although a de- 
signer will not know how close his/her design is to the opti- 
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Figure 2. 
LTJG Ugur Ozkan Using HAMLET. 








mum, our study indicates that deviations are probably less than 
10% away. Ours was the first systematic study of the problem. 

Our experience from this analysis showed the need for a 
tool to analyze PLA design heuristics. in 1988, we launched 
an extensive design project for such a tool. The result was 
HAMLET (Heuristic Analyzer for Multiple-Valued Logic Ex- 
pression Translation), the first practical computer-aided 
(CAD) tool for multiple-valued logic. HAMLET"® has the 
ability to analyze heuristics specified by the user, using either 
individual functions or randomly generated functions. This is 
a particularly useful feature. Instead of writing a separate 
program to generate test inputs, this is done automatically in 
HAMLET. Of special significance is HAMLET’s ability to 
collect and automatically display heuristics. Tasks that nor- 
mally took man-days can now be literally done in keystrokes. 
HAMLET can produce ensembles of test cases whose statis- 
tical properties can be controlled. The search algorithm used 
in reference 13 was also incorporated into HAMLET. We are 
now able to use the search to find solutions that were impossi- 
ble to obtain with any previous method. The search can be 
controlled. At one extreme, it can go directly to a PLA realiza- 
tion of a given function. At the other, itcan perform a complete 
search, producing a known minimal solution. In effect, this is 





like a dial, one end of which corresponds to fast but poor 
solutions, while the other end corresponds to slow but minimal 
solutions. HAMLET can also do design, producing the layout 
of a PLA, given a specification. Figure 3 shows one PLA 
layout produced by HAMLET. HAMLET is designed to have 
a long lifetime. It is easily modified; indeed an improved 
heuristic has been added‘, as well as an entirely new ap- 
proach* (see below). Also, HAMLET has been used by re- 
searchers in at least four countries. 


Simulated Annealing 


When metal or glass is heated and then slowly cooled, it 
is said to undergo annealing. The final cooled state is similar 
to a crystal and represents a low energy state. If it is cooled at 
a high rate, atoms have less of a chance to align and a higher 
energy state occurs. 

In 1953, Metropolis'’ proposed a method of optimizing 
complex problems that mimics this process. It is called simu- 
lated annealing. We have formulated simulated annealing for 
the PLA minimization problem as follows. One can satisfy a 
given PLA specification by producing a set of implicants. Each 
implicant represents a part of the whole PLA. Solving the PLA 
minimization problem, therefore, is equivalent to finding the 
fewest implicants that satisfies the specification. Indeed, there 
is a direct relationship between the number of implicants and 
the chip area occupied. Thus, reducing the number of im- 
plicants by half, for example, results in almost a one-half 
reduction in chip area. Implicants can be manipulated. For 
example, one implicant can be divided into two, pairs some- 
times can be combined into a single implicant, and pairs can 
sometimes be changed into another pair, or a triple, etc.. 
Ideally, one would like to combine a pair of implicants into 
one, since this reduces PLA chip area. However, in a given set 
of implicants, it is quite possible that no pairs combine. Such 
a set is said to be a local minimum. However, by replacing a 
pair by another pair, or by a triple, etc., it is possible to create 
a new set of implicants, where now many combine, leading to 
a smaller set. The descriptor “local” means there are no im- 
plicant sets nearby that contain fewer elements. It is not a 
“global” minima if there is another set of implicants (some- 
where) that has fewer elements. 

The derivation of a minimal set is a difficult problem. A 
heuristic generates a solution to a problem, but not necessarily 
a minimal one. The PLA minimization problem belongs to the 
set of NP-complete problems, which means that the best 
known algorithm for finding the minimal solution increases 
exponentially in complexity as the problem size increases. It 
is for this reason that there is so much interest in heuristics. 
Simulated annealing has the benefit that, under certain condi- 
tions, a minimal solution will be found with a probability that 
approaches 100% as the process continues. This is unlike any 
other known heuristic for PLA minimization, where, for cer- 





Figure 3. 


PLA Layout Produced by HAMLET. 








tain PLA specifications, the probability of finding a minimal 
solution is 0%. 

In effect, simulated annealing is a search through a space 
of solutions. The basic computation is a move. That is, given 
one solution to the PLA specification, another solution is 
obtained by making a move. Each move begins by first select- 
ing a pair of implicants. If the two implicants can be combined, 
they are, completing the move. If the implicants cannot be 
combined, a replacement is considered with the fewest number 
of equivalent implicants. In effect, a weighted coin is tossed. 
If the coin lands with heads up, the replacement is made, 
completing the move. Otherwise, it is not made, and another 
pair is selected. The analogy with annealing is this. At each 
temperature, the weighting of the coin is the same for some 
fixed number of moves. In our programs, we allow 20,000 
moves at each temperature. At high temperatures, the weight 
is such that heads almost always comes up, which means that 
all moves are accepted. However, at low temperatures, the 
weight is shifted to tails, so that fewer moves are accepted 
which increase the number of implicants. At any temperature, 
a move that reduces the number of implicants is accepted. 
Thus, at low temperatures, the mixture of generated moves is 
predominantly PLA complexity-decreasing moves. The exis- 
tence of complexity-increasing moves at all temperatures is 
important; it allows the system to escape from local minima. 
This is analogous to annealing, in which atoms at high tem- 
perature move about freely. At lower temperatures the freedom 
is restricted, as they settle into a low energy state. 

Figure 4 shows how simulated annealing performs on a 
multiple-valued function that begins with 96 implicants. This 
solution was obtained by applying a heuristic to some given 
PLA specification. The axis labeled “number of product 
terms” measures the quality of the solution to the PLA speci- 
fication; the fewer the better. The axis marked temperature 
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Figure 4. 


Simulated Annealing Applied to Multiple-Valued PLA Minimization. 
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measures the probability with which implicant increasing 
moves are accepted. At higher temperatures (toward the front), 
almost all implicant-increasing moves are accepted. The ver- 
tical axis represents the number of times a solution is found at 
the temperature and the number of implicants on the x-y plane. 
Here, color encodes the vertical deviation. In effect, the verti- 
cal axis shows how the simulated annealing program is doing 
with respect to the number of product terms. Specifically, at 
the highest temperature (the slice closest to the front), one can 
see how the simulated annealing progresses from the initial 
solution of 96 implicants up to as many as 216 implicants. This 
is described as melting. As the temperature decreases from this 
point, there is steady progress towards PLA specifications with 
fewer implicants. In this application of simulated annealing, a 
solution of 88 implicants is achieved. The behavior of the 
implicants is similar to the behavior of individual atoms as 
metal or glass anneals; as the temperature decreases, the form 
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of the final material takes shape as individual atoms occupy a 
site in the structure corresponding to low energy. 

A comparison of simulated annealing with other PLA 
minimization heuristics shows that it is superior on an average 
case basis. That is, using HAMLET, random functions were 
generated, heuristics were applied to all, and the results com- 
pared. Simulated annealing surpasses the all known heuristics 
on the solutions obtained. Although, more time is needed, 
simulated annealing has the benefit that its performance can 
be controlled. That is, with larger computation times, one can 
achieve a more compact realization. 


Conclusions 


The Navy has a tradition of commitment to computing, as 
epitomized by the contributions of Rear Adm. Grace Murray 
Hopper. Indeed, in the past 50 years, significant improvement 





has occurred in Navy systems because of the computer. It is 
Clear that the next 50 years will bring even further improve- 
ment. The promise of compact, high speed computers has 
inspired our work on multiple-valued logic. 

Multiple-valued logic offers a means to overcome signif- 
icant disadvantages of binary. It makes more efficient use of 
chip area by encoding more information in the wires that 
connect devices, and it allows the use of high-speed carry-less 
operations. Our work at the Naval Postgraduate School has 
been focused on development of techniques (PLA’s) and the 
tools (CAD) needed to make multiple-valued logic a reality. 

Our current effort focuses on expanding simulated an- 
nealing. This includes a promising approach to high-speed 
computations; indeed, initial indications show a modification 
of the conventional simulated annealing paradigm can pro- 
duce computation speeds higher than any known heuristic, 
while maintaining almost the same capability to minimize chip 
area. We are also pursuing the use of parallel computers to 
design PLA’s®, including simulated annealing. 
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Subsonic and Supersonic 
Mixing and Combustion 
Enhancement 


Gabriel D. Roy 
Office of Naval Research 


Abstract 


A five year research program was sponsored by the US 
Department of the Navy to obtain the science base for the 
conceptual evaluation of higher speed, long range supersonic 
combustion ramjets. This program was focused to determine 
the turbulent structure, size scales, intensities and rates of heat 
and mass transfer of compressible subsonic and supersonic 
shear layer mixing flows, with and without combustion. Al- 
though the program addressed the issues relating to the com- 
bustion of solid fuels in supersonic combustion as well, the 
present paper will focus on supersonic shear layer mixing and 
combustion enhancement. The studies showed that heat re- 
lease tends to inhibit the mixing of shear layers. However, 
three-dimensional excitation at subharmonics of the shear 
layer, and streamwise vorticity stirring by lobbed splitter 
plates have been found to increase mixing and combustion. 
The supersonic shear layers are found to be stable, and their 
growth deceases, as the convective Mach number (Mc-) in- 
creases. Vortex generators increase the growth rate of super- 
sonic shear layers by about 30%, whereas shock impingement 
has very little effect. Non-axisymmetric nozzles and inlets 
have shown to enhance mixing and combustion of supersonic 
coaxial flows. 


Introduction 


As longer ranges are desired for air launched tactical 
missiles, renewed interest has emerged on air-breathing pro- 
pulsion systems. Aircraft and ship-imposed constraints on 
length, diameter and weight of weapons require compact 
ramjets and scramjets. However, scramjet combustors have 
been suffering decreased mixing associated with supersonic 
compressible shear layers, and the associated difficulty in 
achieving complete combustion in the limited residence times 
available. To address the fundamental supersonic mixing and 
combustion processes, the Office of Naval Research (ONR) 
initiated a focused five year research program. The experimen- 
tal, analytical and numerical research program addressed the 
conceptual evaluation of higher speed, longer range super- 
sonic combustion ramjets and the potential for several fold 
enhancement of the radiant output of high performance air- 
craft-launched pyrotechnic infrared countermeasures. In par- 
ticular, the initiative investigated the temporal and spatial 
scales, and rates of subsonic and supersonic mixing and com- 
bustion processes, the supersonic flame structures in mixing 
layers adjacent to bounding solid surfaces, and the character- 
istics of the gasification and particle ejection processes from 
high energy polymeric solid fuels and their boron/metal filled 
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composites into the supersonic mixing layer and flame struc- 
ture. Among others, the program addressed the following 
issues, which are specific to this paper: 








1. | What are the turbulent intensities and size scales in 
compressible subsonic and supersonic shear layer mix- 
ing flows? Are the structures large scale and coherent? 


Do the streamwise narrow band velocity fluctuations, 
momentum thickness and mixedness of these turbulent 
shear layers behave similarly with streamwise distance 





Figure 2. 


Temperature Distribution in the Mixing Layer with and without 
Energy Release 




















as in those compressible subsonic shear layer mixing 
flows? 


Can the supersonic mixing layers be manipulated by low 
level forcing at subharmonics of the natural frequency 
of the mixing layers? Are there other methods of con- 
trol? 

Can the compressible subsonic and supersonic mixing 
layers and flow structures be influenced by combustion 
generated heat release? 


This paper is limited to the research conducted under this 
special focus program, and complimentary research under the 
ONR core program. The approach undertaken to solve the 
issues, and the results obtained by the investigators are sum- 
marized. The extensive detailed information obtained from 
this program is omitted on purpose. The reader may refer to 
the references given at the end of the paper for further infor- 
mation. 


Effect of Heat Release 


To understand the dynamics and topology of the coherent 
structures controlling the development of reactive shear flows, 
in order to develop both a conceptual model and an analytic 
framework for their description, Grinstein performed numer- 
ical experiments. His approach is to use direct and large-eddy 
numerical simulation (LES) of spatially evolving three-di- 
mensional compressible shear flows with chemical reaction 
using monotone 3D flux-corrected transport (FCT) algo- 
rithms’. The details of the Reactive Shear Flow Model and 
the incorporation of chemistry can be found in Ref. 3. The 





Figure 3. 


Mean Profiles of Product Concentration and Temperature at 
Two Streamwise Locations. 
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time-dependent, compressible, reactive flow conservation 
equations for total mass and energy density, momentum, and 
chemical species number densities are solved. The term asso- 
ciated with molecular diffusion, viscous diffusion and thermal 
conduction in the energy equation are calculated explicitly 
using central-difference approximations and coupled to con- 
vection using timestep-splitting techniques. Since solving a 
detailed set of chemical kinetic rate equations in conjunction 
with the unsteady, three-dimensional flow equations is beyond 
current computer resources, simplified global-chemistry mod- 
els were used to understand the effects of chemical energy 
release on the complex three-dimensional flow field. 


Reactive Shear Flow Simulations 


Hussain and Hussain showed that coherent structures are 
capable of enhancing or reducing the mixing of coflowing 
reacting streams, and hence can lead to methods of controlling 
combustion efficiency through control of the initiation, evolu- 
tion, and interaction of coherent structures*. Chemical reac- 
tions between the coflowing streams may also affect the 
structure of the mixing layer, either through compressibility 
effects such as baroclinic vorticity production and expansion 
as energy is released, or by changing the frequency content of 
the mixing layers with the introduction of characteristic fre- 
quencies of the chemical processes. The combustion time 
scales themselves, in turn, may be modified by the natural 
acoustic frequencies of the mixing layer or by those imposed 
through excitation. The interaction between these effects can 
have dramatic consequences on combustion efficiencies. Ear- 
lier studies™* have indicted that the energy released by the 
chemical reaction has the effect of reducing the entrainment 
in the mixing layer. 

The results of numerical simulations of spatially evolving, 
chemically reacting mixing layers are presented. The system 
studied consists of a mixing layer formed by dilute non- 
premixed subsonic streams of hydrogen and oxygen, which 
are allowed to react both with and without energy release. 
Instantaneous results from reactive two-dimensional mixing 
layer simulations are shown in Fig. 1. Here, the product 
concentration is plotted against the length of the mixing layer 
without energy release (Figure 1a), and with energy release 
(Figure 1b). The reaction considered is 2H2 + O2 -> 2H20, 
initial temperature = 1400 K, Mach number = 0.3, reactant 
(molar) concentration is 10% , and the velocity ratio is 3. The 
flow was organized by forcing the cross stream velocity at the 
inflow. 

Spanwise vortices roll up with the forcing frequency as a 
result of the nonlinear evolution of the Kelvin-Helmholtz 
instability, and vortex pairings are inhibited. Thus the growth 
of the mixing layer is significant only in the vortex-roll-up 
region. The center of the shear layer gradually shifts towards 
the upper, slower side as the flow moves downstream, indicat- 
ing the inherently asymmetric entrainment in the planar shear 





Figure 4. 


3-D Simulation of Product Generation with Energy Release and 
Spanwise Excitation. 








3D. with energy release 
& spanwise excitation 


“ORGY 


tite 


HHH 

















Figure 5. 
Enhancement of Vorticity and Chemical Reaction by Second- 
ary Instabilities 
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flow’*, Due to the fast flow regimes considered, diffusive 
mixing takes place in the immediate neighborhood of the 
interface between reactants, and the chemical product gener- 
ation is observed mainly away from the cores of the structures, 
especially in its outer edges, which are the most chemically 
reactive regions. As can be seen from Figure 1, chemical 
energy release has the effect of reducing the shear layer 
growth, and the amount of product formed. The product peak- 
value is found to be 30% larger without energy release. The 
temperature distribution with and without energy release is 
shown in Figure 2. It is noted that the chemical energy release 
is accompanied by significant changes in the temperature, with 
peak-temperatures at 1.23T, and 1.02T,. with and without 
energy release respectively. 

In Figure 3, profiles of time-averaged product concentra- 
tion and temperature are shown at a streamwise location x; 





Figure 6. 


Schematic of the Shear Layer Facility. 
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Figure 7. 


Velocity Profiles at Various Axial Locations 
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Figure 8. 


Experimental Values of Shear Layer Spreading Rate. 
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and another stream wise location x2 downstream for cases with 
and without heat release. These profiles clearly show a shift 
of the flame towards its slower side, which becomes more 
pronounced moving downstream. Part of the asymmetry may 
be attributed to the different mass diffusivities on each side of 
the shear layer. 

However, with spanwise excitation, three-dimensional 
simulations have shown an increase in mixing even with heat 
release. For the same conditions used in the previous simula- 
tion, the spatially evolving product concentrations are shown 
for two consecutive time steps in Figure 4. The three-dimen- 
sional nature of the mixing layer and mixing enhancement is 
evident. Spectral simulations were performed by Met Calfe 
and Hussain to understand the topology and dynamics of large 
scale coherent structures, and their coupling with fine scale 
mixing. These simulations have shown that secondary in- 
stabilities enhance mixing by convoluting and stretching 
flame sheet, and they interact with “rolls” to control large scale 
species transport and reaction rate. Figure 5 shows the vortic- 
ity magnitude and instantaneous product generation with 2-D 
excitation, and with 2-D + 3-D excitation. The three-dimen- 
sional excitation produces intense fine scale mixing, as shown 
by the vorticity magnitude and the instantaneous product 
generation’. 


Papamouschou (Ref. > 
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Supersonic Mixing Studies 


Mixing of Parallel Supersonic Streams 


A study was initiated by Waltrup and Sullins to determine 
the structure and growth of nonreacting shear layers formed 
when a sonic or supersonic two-dimensional jet mixes with a 
two-dimensional supersonic air stream. The primary airstream 
of a typical Mach number of 2 or 3 is separated from a sonic 








Figure 9. 


Schlieren Photograph of the Shear Layer. 
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Figure10. 


Normalized Growth Rates and Peak Turbulence Quantities as 
a Function of Relative Mach Number. 
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airstream by a sharp trailing edge splitter plate. This is fol- 
lowed by a constant area supersonic mixing section (Figure 6). 
Diagnostics included Schlieren photography, and pitot and 
cone static pressure and velocity measurements using a one 
component laser Doppler velocimeter. 

If both streams are incompressible (M < 0.3), the shear 
layer is inherently unstable, and large scale, coherent vortical 
structures are formed. The growth is dependent on vortex 
pairing and the ratio of the velocities of the two streams. At 
supersonic speeds, the shear layer tends to be more stable, and 
efficient mixing is not achieved. The concept of convective 
Mach number Mz, established by Roshko and Papamoschou!” 
applies. 


Tests were conducted with two-dimensional shear layers 
produced by both Mach 2 and Mach 3 nozzle streams mixing 
with a Mach 1 nozzle stream. Instream cone static and pitot 
probes were used to determine the velocity profiles at various 
axial locations (Figure 7). The velocity profiles were used to 
determine the shear layer spreading rate (Figure 8). Schlieren 
photographs were also taken to determine the spreading rate 
(Figure 9), and the spreading rate determined from the velocity 
profiles compared well with these photographs. 

In a series of experiments, Krier and his associates inves- 
tigated the fundamental fluid dynamic and combustion mech- 
anisms occurring in the entrainment, mixing, and combustion 
processes in the free shear layer formed between two high 
speed (supersonic/supersonic or supersonic/subsonic) 


Mean and Turbulent Velocity Measurements: Velocity 
measurements have been obtained using a two-component 
laser Doppler velocimeter system in a two stream, supersonic 
mixing layer facility. A variety of cases has been examined 
with free stream velocity ratios in the range from 0.16 to 0.78, 
freestream density ratios ranging from 0.57 to 1.55, and rela- 
tive Mach numbers from 0.40 to 1.97. (The relative Mach 
number is used here as a measure of compressibility effects 
and is defined as the freestream velocity difference divided by 
the mean speed of sound, M, = A U/ a). The range of relative 
Mach numbers examined spans the range of significant com- 
pressibility effects as determined by the reduction in normal- 
ized shear layer growth. 

Measurements from the developing regions of these shear 
layers suggest that a local Reynolds number based on 
freestream velocity difference and mixing layer thickness of 
Reg = A UB/ ¥=1x 10° is required for full development of 
the mean turbulent velocity fields. To determine the growth 
rates, least-squares linear regressions were performed on the 
mixing layer thickness data from the fully developed region. 
The resulting growth rates have been normalized by the 
growth rate of corresponding incompressible mixing layer at 





Figure11. 


2-D Scalar Covariance Field for M; = 0.63. (a) Transverse 
Image, (b) Oblique Image 
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the same freestream velocity and density ratios, and these 
normalized growth rates are plotted in Figure 10 with respect 
to the relative Mach number. The results clearly show a 
reduction in the normalized growth rate due to compressibility, 
which is consistent with those observed by Roshko and others. 

Variations in the peak values of the fully developed high- 
speed mixing layer turbulence quantities from their incom- 
pressible counterparts reflect the effects of compressibility on 
turbulence. The peak values of the streamwise turbulence 
intensity, transverse turbulence intensity, and the normalized 
kinematic Reynolds shear stress, <u’ v’>/ (AU), for the com- 
pressible cases have been normalized by the respective incom- 
pressible values of 0.18, 0.14, and 0.013 and are plotted in 
Figure 10. From order-of-magnitude estimates using boundary 
layer approximations to the governing equations, it can be 
shown that for a compressible mixing layer the normalized 
kinematic Reynolds shear stress should scale with the normal- 





Figure1 2. 


2-D Scalar Covariance Field for M, = 1.49. (a) Transverse 
Image, (b) Oblique Image. 
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Figure13. 


Comparison of Product Generated with and without 3-D Exci- 
tation 
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ized growth rate. The comparable reduction in these two 
quantities with increasing relative Mach number is seen in 
Figure 10. 

Large Scale Structures: The large scale structures of com- 
pressible mixing layers have been studied by Krier’s group 
using planar Mie scattering imaging and statistical analysis of 
two-dimensional covariance fields'*. Three flow conditions 
were investigated, using both traditional transverse imaging 
as well as oblique imaging. Two types of seeding methods 
were employed, passive scaler transport of sub-micron ethanol 
droplets from one freestream to the other as well as the 
formation of ethanol droplets due to molecular mixing of air 
from one stream, laden with ethanol vapor, with cold air from 
the other stream. The relative Mach numbers of the three flow 
conditions were M, = 0.63, 1.24 and 1.49, while the local 
Reynolds numbers for these cases were Re = A Ub/ v = 3.0 
x 10° and 6.3 x 10° respectively. The level of three-dimension- 
ality in a compressible mixing layer was found to be largely a 
function of relative Mach number with the Reynolds number 
apparently having little effect other than to introduce a wider 
range of turbulent structure scales. 

The instantaneous images provide some striking exam- 
ples of the various flow structures within compressible mixing 
layers. Organized large scale structures were present in all 
cases; however they appeared most frequently in the low 
relative Mach number flow. As the relative Mach number is 
increased, the transverse images suggest the structures become 
less organized, with an increased presence of randomly ori- 
ented, smaller scale. Structures indicative of streamwise 
counter-rotating vorticity, and a measure of the level of three-di- 
mensionality of the mixing layer, were more often observed at 
higher relative Mach number cases. The sizes of these structures 
also tended to increase with the frequency of occurrence. 





Figure14. 


Typical Convoluted Splitter Surface Showing Periodic Stream- 
wise Vortex Array. 
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The spatial covariance fields of the ensembles of 256 
instantaneous images for each case were computed to deter- 
mine, on a sound statistical basis, the characteristics of the 
large scale structures, Figure 11 and 12 show the covariance 
fields for both the transverse and oblique scalar transport 
images of the low and high relative Mach number cases, 
respectively. In Figure 11 and 12, the central peak is 100%, 
and each successive contour is 20% lower. These results 
indicate that, as the relative Mach number is increased, the 
transverse scalar transport structures tend to increase in 
dimensionless size ratio and eccentricity while the angular 
orientation decreases slightly. These results agree extremely 
well with recent time-evolving numerical simulations!* of 
compressible mixing layers which also show more flattened, 
elliptical structures at large relative Mach number. Oblique 
scalar transport correlations indicate structures of generally 
smaller size ratio and significantly reduced eccentricity than 
their transversely imaged counterparts, yet the effect of in- 
creased relative Mach number indicates only a slight increase 
in dimensionless structure size with no apparent trend in the 
eccentricity. The product formation studies exhibited struc- 
tures, in both the transverse and oblique image planes, which 
were smaller in dimensionless size and less elliptical in shape 
than the top stream scalar transport structures. Significantly 
different characteristics between transverse and oblique real- 
izations of the scalar transport structure implies the presence 
of an organized, large scale mixing layer structure, which is 
not entirely due to increased three-dimensionality of the mix- 
ing layer associated with increased relative Mach number. 





Figure15. 


Comparison of Flame Spreading Achieved with a Flat Splitter 
and a Convoluted Splitter. 
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Figure16. 


Vortex Ring Penetration with Transverse Pulse Injection 
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Figure17. 


Pitot Tube Growth with Planar Shock Impingement on bound- 
ary Layer. 








—a— DSLI0-BL1 
a-=@-- DSL10-BL2 
DSL20-BL 

















Four!1992 —One/1993 15 





Mixing and Combustion 
Enhancement 


Subsonic Flows 


As mentioned in section II, three-dimensional excitation 
has been found to increase mixing and combustion. The en- 
hancement in vorticity and product generation by excitation, 
shown in Figure 5, is further quantified in Figure 13. Here the 
ratio of product generation with 3-D excitation to that with 2-D 
excitation is plotted against time. As can be seen, over 30% 
increase in product generation is obtained by 3-D excitation. 

Streamwise Vorticity Stirring: Experimental work by 
McVey focussed on the effect of streamwise vorticity on the 
spreading of flames in high speed flows. Using a convoluted 
splitter surface as shown in Figure 14, and by the proper choice 
of the amplitude, pitch and profile of the lobed surface, it was 
possible to generate arrays of very-large-scale, intense vorti- 
ces. Flow visualizations indicated that the influence of these 
vortices is first, to cause a large scale exchange of fluid 
(stirring) followed by a rapid breakdown of the organized 
structures into random small scale turbulence, which is essen- 
tial to high speed combustion "4. 

Non-premixed flame with air as the oxidant and a mixture 
of nitrogen and vaporized triethylborane as fuel was used. 
Since the fuel is pyrophoric, no ignition source was required. 
Flame photographs with a flat splitter and a lobed splitter are 
shown in Figure 15. A dramatic increase in the rate of flame 
spreading is evident, which is attributed to the streamwise 
vorticity stirring. Further work to understand the details of the 
mechanisms involved is presently sponsored by the Army 
Research Office. 





Figure18. 


Pitot Thickness Growth with 3-D Shock Impingement on Shear 
Layer. 
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Figure19. 


Vortex Generators 
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Figure20. 


Pitot Thickness Growth with Cylindrical Vortex Generators. 
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Transverse Pulsed Fuel Injection: Vakili and Wu are studying 
transverse pulse injection of fuel into the main oxidant stream, as 
a means of enhancing mixing and combustion. Experiments 
conducted in subsonic flows indicated that significant penetration 
of the vortex ring can be obtained by transverse pulse injection 
(Figure 16). In particular, sharp edge orifices have been found to 
be very effective. Further studies will be undertaken to determine 
the effectiveness in supersonic mixing enhancement, and conse- 
quently supersonic combustion’®. 


Supersonic Flows 


Compressible shear layers are very stable, and the spread- 
ing rate decreases as the convective Mach number increases. 
In order to achieve complete combustion within the short 
residence times available in supersonic combustors, it is nec- 
essary to enhance mixing. Energy release could further com- 
plicate and decrease mixing. Experimental studies have been 
undertaken by Dolling, and by Schadow to enhance supersonic 
mixing, and mixing and combustion respectively, by passive 
or active means. 

Dolling’s research focused on the growth rate and struc- 
ture of high Reynold’s number turbulent shear layer formed 
by Mach 5 and Mach 3 airstreams'®. The first phase of the 
study was to determine whether the growth rate of the shear 
layer could be increased by changing the initial conditions at 
the shear layer origin. The initial conditions were changed 
either by a planar 3-D shock impingement on the turbulent 
boundary layer upstream of the origin, or by using cylindrical 
or wedge shaped vortex generators in the boundary layer 
upstream of the origin. 

Shock Impingement: Both planar and 3-D shock waves 
were generated and arranged to impinge on either the bound- 
ary layer upstream of the shear layer origin or at the shear layer 
itself. Seven cases were studied. Typical results will be pre- 
sented here. 

For planar shock impingement on boundary layer, three 
different configurations were used. In Figure 17, the growth 





Figure21. 


Experimental Set-up (Co-axial jets) 
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rate of the shear layers for these cases are shown. The growth 
rate of the undisturbed shear layer is also shown, by hatch line, 
for comparison. A maximum of 30% increase in shear layer 
growth rate was obtained. However, within a few boundary 
layer thicknesses downstream, the growth rates are not much 
different from the undisturbed growth rate. 

Three-dimensional shock impingement on the shear layer 
generated only slight three-dimensionality in the layer im- 
mediately downstream of the impingement location due to the 
pressure gradient caused by the sine wave shape of the incident 
shock. However, the flow became more two-dimensional fur- 
ther downstream, and the thickness of the shear layer remained 
about the same. The growth rates across the span are summa- 
rized in Figure 18. As can be seen, there are no observable 
differences between the growth rate and thickness for this case 
with the undisturbed case. 

Vortex Generation: Two types of vortex generators were 
used. The first consisted of a row of small cylinders screwed 
into the test surface, and the second consisted of corrugated 
plates or wedges (Figure 19). The results show that the shear 
layer growth rate can be increased by about 30% using cylin- 
drical vortex generator arrangements (Figure 20), whereas 
with wedge-shaped arrangements, the growth rate is unaf- 
fected or even reduced. 

Non axisymmetric Nozzle Flows: Schadow and his asso- 
ciates have demonstrated the superiority of non axisymmetric 
nozzles, with respect to circular nozzles, in terms of fine scale 
and large scale mixing in free-jet and shrouded-jet tests with 
and without reaction'”'*, They have extended this work to 
supersonic jets. The experiments were performed in a coaxial 
supersonic jet facility, which can be used for nonreacting 
mixing studies and for supersonic fuel-rich plume combustion 
studies (Figure 21). 

For the mixing studies, the inner (center) jet was kept at 
M; = 3. Helium, argon, nitrogen, and heated and unheated air 
were used to change the density of the inner jet. The velocity 





Figure22. 


Normalized Spreading Rate of Circular and Noncircular Jets 
as a Function of Convective Mach Number. 
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Figure23. 


Radiance Contours of Circular and Multi-step Nozzles 
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Figure24. 


Relative Total Radiance for Circular and Multi-step Nozzles 
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of the outer jet was varied from zero (no coaxial flow) to fully 
expanded condition (M2 = 1.8) for ambient and 475K air 
temperatures. Six test conditions were selected to obtain a 
convective Mach Number range of 0.5 M, 2.2. Mixing be- 
tween the center jet and outer air was determined with a 
total-pressure probe, and gas sampling probe connected to an 
analyzer. 

In Figure 22, the variation of the normalized spreading 
rates of the circular jet and the rectangular jets (two configu- 
rations) are shown against the convective Mach number. The 
behavior observed for the circular coaxial jet is similar to that 
of the planar shear layer, with the normalized spreading rate 
decreasing with increasing convective Mach number. How- 
ever, the spreading rates of the rectangular jets are higher than 
the circular jet, both in the major and minor axis planes, in 
almost the entire convective Mach number range. The spread- 


18 Naval Research Reviews 


ing rates are relatively higher for M. 0.4, because the rectan- 
gular jets have a higher incompressible spreading rate than the 
circular jet. 

The supersonic fuel-rich plume combustion experiments 
were also conducted in the same coaxial jet facility. Here, the 
outer flow at ambient temperature was fixed at M; = 1.3 and 
the inner fuel-rich flow at M; = 2. The fuel-rich plume was 
produced in a gas generator by burning ethylene, air, and 
oxygen at an equivalence ratio of 2, and a combustion temper- 
ature of T; = 2300K. Combustion characteristics are compared 
based on the infrared images, which indicate the spreading rate 
and combustion intensity. As a typical example, in Figure 23, 
the radiance map of the underexpanded coaxial jets from a 
circular and multi-step nozzle are compared '*. The RMS level 
of the multi-step nozzle is higher than the circular jet. In Figure 
24, the circular and multi-step nozzles are compared for three 
pressure ratios based on the relative integrated radiance from 
the nozzle lip to x/D- = 16. The radiance is normalized by the 
value for the circular jet at the design pressure ratio. These 
figures indicate the superior performance of the multi-step 
nozzle, and the increased heat release associated with en- 
hanced fine-scale mixing. 


Summary 


A special focus program was sponsored by ONR to un- 
derstand the mechanisms involved and the control methodol- 
ogy to enhance mixing and combustion is subsonic and 
supersonic reacting and nonreacting flows. This multi-disci- 
plinary, numerical and experimental research by the academia, 
Navy and industry laboratories produced the following signif- 
icant results: 


1. Energy release has a tendency to reduce mixing and 
further product generation. 


Three-dimensional excitation increases product genera- 
tion (with energy release). 


Secondary instabilities are very important in both the 
core and braid regions of the flow, and they interact with 
the rolls to control large scale species transport, and thus 
the reaction rate. 


Streamwise vorticity stirring by means of convoluted 
splitter plates has shown a significant improvement in 
flame spreading. 


Transverse pulsed injection of the fuel allows deep 
penetration of the vortex rings into the mainstream, and 
enhances mixing. 


Supersonic mixing layers are very stable, and the growth 
of the mixing layer decreases as the convective Mach 
number increases. 


Shock impingement in the boundary layer in supersonic 
flow produces a small increase in the shear layer growth, 





but within a few boundary layer thicknesses down- 
stream, it becomes similar to the undisturbed growth. 


Supersonic shear layer growth is increased by about 
30% by vortex generation, using cylindrical vortex gen- 
erators. 


Non axisymmetric nozzles have increased mixing rate 
as compared to circular nozzles. In coaxial supersonic 
flows, they have produced increased fine scale mixing, 
and enhanced combustion. 


Further research and development is needed to control 
and enhance supersonic mixing and combustion. Active, pas- 
sive, and hybrid control methodologies may be required for 
specific application. ONR is currently pursuing another spe- 
cial focus program on active control of complex flows. 
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Models for Assessing the 
Reliability of Computer 


Software 


Nozer D. Singpurwalla and Simon P. Wilson 
The George Washington University 


Abstract 


A formal approach for evaluating the reliability of com- 
puter software is through probabilistic models and statistical 
methods. This paper is an expository overview of the literature 
on the former. The various probability models for software 
failure can be classified into two groups; the merits of these 
groups are discussed and an example of their use in decision 
problems is given in some detail. The direction of current and 
future research is contemplated. This research has been sup- 
ported by the Office of Naval Research, the Army Research 
Office, and the Air Force Office of Scientific Research. 


1. Introduction 


Having been developed over the last 20 years, software 
reliability is a relatively new area of research for the statistics 
and the computer science communities. It arose because of 
interest in trying to predict the reliability of software, particu- 
larly when its failure could be catastrophic. Obviously the 
software that controls an aircraft carrier, a nuclear power 
Station, a submarine or a life-support machine needs to be very 
reliable, and statistical techniques will aid the computer scien- 
list in deciding if such software has sufficient reliability. The 
subject is also of commercial importance, as for example when 
decisions have to be made concerning the release of software 
into the marketplace. 


All software is subject to failure, due to the inevitable 
presence of errors (or bugs) in the code, so the first aim of the 
subject has been to develop models that describe software 
failure. There are various methods of specifying such failure 
models, and Section 2 discusses these in some detail. It is fair 
to say that this model derivation has been the focus of research 
so far. Once a failure model has been specified then it can be 
applied to problems such as the optimal time to debug software 
or deciding whether software is ready for release. These appli- 
cations have received less attention in the literature but are 
becoming more prevalent. We will mention here that there is 
another approach to software reliability that differs consider- 
ably from the statistical ideas presented here. This approach 
attempts to prove the reliability, or correctness, of software by 
formal means of proof, just as one would prove a mathematical 
theorem. This is an exercise in logic, albeit a rather complex 
one. It works well on small programs, for example on a 
program that computes the factorial function, but becomes a 
forbidding task for even moderately complex pieces of code. 
Nevertheless, the idea that software can be proved correct is 
appealing. The approach is not discussed further. 

This paper is divided into 5 further sections. Section 2 
categorizes the different strategies that have been used to 
model software failure. Section 3 reviews the historical devel- 
opment of the subject by describing some of the more com- 
monly used models, and Section 4 shows that many of these 
models can be unified if one adopts a Bayesian position. 
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Figure 1. 


Categorization of Software Reliability Models. 








Section 5 looks at applications of the material developed in 
Section 3, and Section 6 concludes with a look at the current 
and future direction of the subject. We assume that the reader 
has some familiarity with some basic reliability and probabil- 
ity concepts; in particular it is important that he or she has 
knowledge of some common probability distributions, statis- 
tical inference and decision making, Poisson processes and the 
concept of a failure rate. 


2. Model Categorization 


All statistical software reliability models are probabilistic 
in nature. They attempt to specify the probability of software 
failure in some manner. In looking through the literature, one 
observes that the models developed so far can be broadly 
classified into two categories. 


Type I: Those which propose a probability model for 
times between successive failure of the software, and 

Type Il: Those which propose a probability model for 
the number of failures up to a certain time. Time is often taken 
to be CPU time, or the amount of time that the software is 
actually running, as opposed to real time. In theory, specifica- 
tion via one of these two methods enables one to specify the 
other. So a model that specifies time between failure will also 
be able to tell you the number of failures in a given time, and 
vice versa. In practice, this may not be straightforward. 

The first of these categories, modeling time between 
failure, is most commonly accomplished via a specification of 
the failure rate of the software as it is running. When this is 
the case the model is to be of Type I-1. The failure rate for the 
i-th time between failure is given, for i=1, 2, 3,... and a 
probability model results. One distinctive feature of software 
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is that its failure rate may decrease with time, as more bugs 
are discovered and corrected. This contrasts with most me- 
chanical systems which will age over time and so have an 
increasing failure rate. An attempt to debug software may 
introduce more bugs into it, thus tending to increase the failure 
rate, so the decreasing failure rate assumption is somewhat 
idealized. However, most of the models of this type that are 
reviewed here have a decreasing failure rate. 

Another way to model time between failure is to define a 
stochastic relationship between successive failure times. Mod- 
els that are specified by this method are known as Type I-2, 
and have the advantage over Type I-1 in that they model the 
times between failure directly, and not via the abstract concept 
of a failure rate. For example, let T,, T.,, ..., T,, ... be the length 
of time between successive failure of the software. Asa simple 
case, one could declare that T,,, = pT, + €,, where p 2 Oisa 
constant and €, is an error term (typically some random vari- 
able with mean 0). Then p < 1 would indicate decreasing time 
between failure (software reliability expected to become 
worse), p = 1 would indicate no expected change in software 
reliability whilst p > 1 indicates increasing times between 
failure (software reliability expected to improve). Those fa- 
miliar with time series will recognize the relationship in this 
example as an auto-regressive process of order 1; in general, 
one would say Thy = {(T,, T,, .. T)) + €; for some function f. 

The second of these categories, modeling the number of 
failures, uses a point process to count the failures. Let M(t) be 
the number of failures in the software that are observed during 
time [O,t). M(t) is modeled by a Poisson process, which is a 
stochastic process with the following properties. 


i) M(0) = 0 and if s <t then M(s) < M(t). M(t) takes 
values in {0, 1,2,...} 





The number of failures that occur in disjoint time inter- 
vals are independent. So, for example, the number of 
failures in the first 5 hours of use has no effect on the 
number of failures in the next 5 hours. 


The number of failures to time t is a Poisson random 
variable with mean j1(t), for some non-decreasing func- 
tion ,(t); that is to say: 


P(M() =n) = HOF --no n=0,1,2.... 


The different models of this type have a different function 
u(t), which is called the mean value function. The mean number 
of failures at time t is indeed ,1(t), as is the variance. The Poisson 
process is chosen because in many ways it is the simplest point 
process yet it is flexible and has many useful properties that can 
be exploited. This second approach has become increasingly 





Figure 2{(i). 


The failure rate of the model of Jelinski and Moranda 
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Figure 2{(ii). 


The failure rate of the model of Littlewood and Verall. 
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popular in recent years. M(t) can also be specified by its 
intensity function A(t), which is the derivative of (t) with 
respect to t; either of these functions completely specify a partic- 
ular Poisson process. One disadvantage of this approach is that it 
implies that there are conceptually an infinite number of bugs in 
the program, which is obviously impossible for code of a finite 
length. Another disadvantage is more subtle; the model implies a 
positive correlation between the number of failures in adjoining 
time intervals, a situation which is not true since again the total 
number of bugs has to be finite. 

Figure 1 is a flow-chart showing the above categorization 
of the statistical models. 


3. Review of Some Software 
Reliability Models 


This section introduces some of the more well known 
probability models for software reliability. There are examples 
of models from each of the two main categories that were 
discussed in the previous section. Since the main purpose of 
the review is to describe the ideas and assumptions behind the 
models, technical details will be kept to a minimum in most 
cases. Those interested in the details of a particular model are 
advised to reference the papers where they were originally 
presented. 

Some common notation will be assumed throughout this 
section and is given below: 


i) Ti= i-th time between failure of the software [i.e. time 
between (i-1)th and i-th failure]. 


ii)  *T,= failure rate for T;, the i-th time between failure, at 
time t. 


iii) M(t)= number of failures of the software in the time 
interval [0, t) (a Poisson process). 





Figure 2{(iii) 


The failure rate of the model of Moranda. 
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A(t)= intensity function of M(t). 


j(t)= expected number of failures of software in time 
[0,t). 


= f A (s) ds, since M(t) is a Poisson process. 
0 


10 models are presented. Model numbers 1 to 7 are of 
Type I-1, models 8 and 9 are of Type II and model 10 is of type 
I-2. Acommon problem to all the models is the lack of data on 
which to test their validity; data on software reliability is 
commercially sensitive and so statisticians in academia have 
very little information on how software in the marketplace 
actually performs. For this reason it is important that the 
assumptions made in deriving these models are carefully 
thought about. 


The Model of Jelinski & Moranda 
(1972). 


This was the very first software reliability model that was 
proposed, and has formed the basis for many models devel- 
oped after. It is a Type I-1 model; it models times between 
failure by considering their failure rates. Jelinski and Moranda 
reasoned as follows. Suppose that the total number of bugs in 
the program is N, and suppose that each time the software fails, 
one bug is corrected. The failure rate of the i-th time between 
failure, T;, is then assumed a constant proportional to N-i+1, 
which is the number of bugs remaining in the program. In 
other words 


fr, (@ IN, A)=A(N-i+ 1), i=1,2,3,...,¢20, 


for some constant A . 


There are some criticisms that one could make of the 
model. It assumes that each error contributes the same amount 
to the failure rate, whereas in reality different bugs will have 
different effects. It also assumes that every time a fix is made, 
no new bugs are introduced; note [see Figure 2(i)] that the 
successive failure rates are indeed decreasing. A model like 
this is sometimes referred to as a "de-eutrophication model", 
because the process of removing bugs from software is akin 
to the removal of pollutants in rivers and lakes. 


Bayesian Reliability Growth Model 


(Littlewood & Verall (1973)). 
Like the Jelinski & Moranda model, the model proposed 


by Litthewood and Verall looked at times between failure of 
the software. However, they did not develop the model by 





Figure 2(v). 


The intensity function for the model of Goel and Okumoto. 
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Figure 2{iv). 


The failure rate of the model of Schick and Wolverton 


Figure 2(vi). 


The intensity function for the model of Musa and Okumoto. 
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characterizing the failure rate; rather they stated that the model 

should not be based on fault content (as Jelinski & Moranda 

had assumed) and then declared that T, has an exponential 

distribution with scale A., and that A, itself has a gamma 

distribution with shape a and scale w(i), for some function wy. 

Despite this it is still considered to be a Type I-1 model. 
Specifically: 


£(tlA)=A.e“% t20 


vy (Wi) )* 40-1. -w@a 
1, Alay (i))= T (a) we A2=0 


w(i) was supposed to describe the quality of the program- 
mer and the programming task. As an example, they chose 
w(i)= B, + B,i. One can show that this makes the failure rate 
of each T; decreasing in tand that each time a bug is discovered 
and fixed there is a downward jump in the successive failure 
rates; see Figure 2(ii). In fact 


a 
——— > 
ry (t1a,B,B,) B+ Bi for t2 0. 


If B, >1 then the jumps in the failure rate decrease in i, if 
B, <1 they increase whilst if B, = 1 they remain a constant. So 
if B, differs from 1 then the fixing of each bug is making a 
different contribution to the reduction in the failure rate of the 
software, an apparent advantage over the model by Jelinski & 
Moranda. This model has received quite a lot of attention and 
has been the subject of various modifications; see models 6 
and 7 later in this section. 


The De-Eutrophication model of 
Moranda (1875). 


Another (de-eutrophication) model of Moranda (1975) 
attempted to answer some of the criticisms of the Jelinski & 
Moranda model, in particular the criticism concerning the 
equal effect that each bug in the code has on the failure rate. 
He hypothesized that the fixing of bugs that cause early 
failures in the system reduces the failure rate more than the 
fixing of bugs that occur later, because these early bugs are 
more likely to be the bigger ones. With this in mind, he 
proposed that the failure rate should remain constant for each 
T;, but that it should be made to decrease geometrically in i 
after each failure i.e. for constants D and k. 


r, (tIDk)=Dk' t 20, D>@ and0<k<1. 


Compared to the Jelinski & Moranda model, where the 
drop in failure rate after each figure was always A, the drop in 
failure rate here after the i-th failure is D k*'(1-k) see Figure 
2(iii). The assumption of a perfect fix, with no introduction of 
new bugs during the fix, is retained. 


Imperfect Debugging Model 
(Goel & Okumoto (1978)). 


This model is another generalization of the Jelinski & 
Moranda model which attempts to address the criticism that a 
perfect fix of a bug does not always occur. Goel & Okumoto’s 
Imperfect Debugging Model is like the Jelinski & Moranda 
model, but assumes that there is a probability p, 0 < p < 1, of 
fixing a bug when it is encountered. This means that after i 
faults have been found, we expect i x p faults to have been 
corrected, instead of i. Thus the failure rate of T, is 


r, (tIN,A,p) = A(N- p (i- 1)) 


When p =1 we get the Jelinski & Moranda model. 


A model by Schick & Wolverton (1978). 


This is yet another Type I model, and this time the failure 
rate is assumed proportional to the number of bugs remaining 
in the system and the time elapsed since the last failure. Thus 


r, (tINA)=A(N-i+1)t, 20 


This model differs from models 1-4 in that the failure rate 
does not decrease monotonically. Immediately after the i-th 
failure, the failure rate drops to 0, and then increases linearly 
with slope (N-i) until the (i+1)th failure; see Figure 2(iv). 


Bayesian Differential Debugging Modei 
(Littlewood (1980)). 


This model can be considered as an elaboration of model 
2 proposed by Littlewood & Verall. Recall that in model 2 it 
was assumed that A,, the failure rate of the i-th time between 
failures, was declared to have a gamma distribution. In this 
new model Littlewood supposed that there were N bugs in the 
system (a return to the bug counting phenomenon), and then 
proposed that A, be specified as a function of the remaining 
bugs. In particular, he stated A. = $, + 9,+...+ $y, where 
, were independent and identically distributed gamma ran- 
dom variables with shape a and scale B. This implied that A, 
would have a gamma distribution with shape o(N-i) and scale 
B. In other respects its assumptions are identical to the original 
Littlewood/Verall model. 


Bayes Empirical Bayes or Hierarchical 
Model (Mazzuchi & Soyer (1988)). 


In 1988 Mazzuchi & Soyer proposed a Bayes Empirical 
Bayes or Hierarchical extension to the Littlewood & Verall 
model. As with the original model, they assumed T, to be 
exponentially distributed with scale A,. Then they proposed 
two ideas for describing A., here called model A and model B. 
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Model A 


Still assume that A. is described by a gamma distribution, 
but with parameters a and B. Now assume that a and B are 
independent and that they themselves are described by prob- 
ability distributions; a by a uniform and B by another gamma. 
In other words: 


n, Q.1o,B)= Zo Att eB A>0 


O<a<v 


mB 0.b) = BY, B>0,a>0,b>0. 


Model B 


Assume that A, is described exactly as in Littlewood and 
Verall i.e. 


vy (WG)? 40-1 --wa@a 
M1, Alay (i) = r (a) zz @ . aes 


and that y(i)= B. + Bri, except now place probability distribu- 
tions on a, Bo and B; as follows: 
1 


nal w)=—., 


0<a<s@ 


2 


m(By | a.b.B,) = l(a) 


(B, + B*- le 6, + BD 
B,2-B,,a>0,b>0 


< 


t (6, c,d) = To 


B, c-1 eB, 


B, 20,c>0,d>0. 


So a is described by a uniform distribution, B, by a shifted 
gamma and B, by another gamma, and there is dependence 
between B, and B,. By assuming B, to be degenerate at 0, 
model A is obtained from model B. The authors were able to 
find an approximation to the expectation of T.,. given that 
T,=t,, T=, ..., T=t,, and so use their model to predict future 
reliability of the software in light of the previous failure times. 


Time-dependent Error Detection Model 
(Goel & Okumoto (1979)). 


This is the first Type II model that we will consider. It 
assumes that M(t), the number of failures of the software in 
time [0,t), is described by a Poisson process with intensity 
function given by 


A(t)=abe™ 
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where a is the total expected number of bugs in the system and 
b is the fault detection rate; see Figure 2(v). Thus the expected 
number of failures to time t is: 


w=] abe ds =a(1-e. 


The function p(t) completely specifies a particular Pois- 
son process, and the distribution of M(t) is given by the well 
known formula 


P(M (t)=n)= f oF en HO, n=0,1,2.... 


Experience has shown that often the rate of faults in 
software increases initially before eventually decreasing, and 
so in Goel (1983) the model was modified to account for this 
by letting 


A(t) =abe fe 


where a is still the total number of bugs and b and c describe 
the quality of testing. 


Logarithmic Poisson Execution Time 
Model (Musa and Okumoto (1984)). 


The Logarithmic Poisson Execution Time Model of Musa 
and Okumoto is one of the more popular software failure 
models of recent years. It is a type II model, but the model is 
not derived by directly assuming some intensity function A(t), 
as was the case with the model of Goel & Okumuto. Here A(t) 
is expressed in terms of 1(t), the expected number of failures 
in time [0,t) via the relationship. 


A(t) =A, ORO. 


Put simply, this relationship encapsulates the belief that 
the intensity (or rate) of failures at time t decreases exponen- 
tially with the number of failures experienced, and so bugs 
fixed, up to time t. The fixing of earlier failures will reduce 
A(t) more than the fixing of later ones. Since we are modeling 
the number of failures by a Poisson process, then we have 
another relationship between A(t) and 1(t), namely 


yt (t) = 4 A (s) ds. 


Using these two relationships between A(t) and 1(t), there 
is a unique solution for the two functions: 


A(t) = M (=F In, O+ 1). 


a I 
A, Ot+ 1 
Figure 2(vi) shows a plot of A(t) versus t; it is similar to 
the plot of figure 2(v) except that the tail is thicker. 
It now follows from the above that by using 
P (M (t) =n) = (u(t) */n! we can say 





(In (A, 6t + 1))" 


P(M = = 1 ’ 
cei 6 (A, Ot+ 1)* xn! 





a=Q(@, 1,2... 


As a final remark, we mention that in their paper the 
authors go into some detail on estimation of A, and @ by 
maximum likelihood methods; however, one of the likelihoods 
appears to be incorrect. 


Random Coefficient Autoregressive 
Process Model (Singpurwalla & Soyer 
(1985)). 


This is a Type I-2 model, that is one that does not consider 
the failure rate of times between failure. Instead it assumes that 
there is some pattern between successive failure times and that 
this pattern can be described by a functional relationship between 
them. The authors declare this relationship to be of the form 

T.=T.% , i=1,2,3,... 
i i-l 
where T, is the time to the first failure and 6; is some unknown 
coefficient. If all the 6;’s are bigger than 1 then we expect 
successive lifelengths to increase, and if all the 6;’s are smaller 
than 1 we expect successive lifelengths to decrease. 


Uncertainty in the above relationship is expressed via an 
error term 6,, so that 


T, =5,T,_,°. 

The authors then make the following assumptions, which 
greatly facilitate the analysis of this model. They assume the 
T,’s to be lognormally distributed, that is to say that log T,’s 
have a normal distribution, and that they are all scaled so that 
T, 21. The 8,’s are also assumed to be lognormal, with median 
1 and variance G,” (the conventional notation is A(1, 6,”)). 


Then by taking logs on the relationship above they obtain 
log T, = 8, log T,_, + log 5. 


=6. logT,_,+€, ,say. 





Figure 3. 


Decision process for single-stage testing. 
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Since the T,’s and the 8,’s are lognormal so the log T.’s 
and the €,’s (= log 5,’s) will be normally distributed, and i in 
particular: €, has mean 0 and some variance o? (the conven- 
tional notation is N(0, o”)). The log-lifelengths therefore form 
what is known as an autoregressive process of order 1 with 
random coefficients 8.. There is an extensive literature on such 
processes which can now be used on this model. 

All that remains to do is to specify 9, and the authors 
consider several alternative models. For example, one could 
make 6. itself an autoregressive process: 

§.=a8.,+@, , where w,is N (O, W.) with W, known. 

When @ is known, the expressions for log T,, and 86. 
together form a Kalman filter model, on which there is also an 
extensive literature. When o is not known the solution is via 
an adaptive Kalman filter algorithm for which the authors 
propose an approach. As an alternative to the above, one could 
place a two stage distribution on @,, and the authors considered 
the idea of 0, being N(A, o, 2). with 2 also a normal random 
variable having mean m, and variance s 2 In this latter case 
one can employ standard hierarchical Bayesian inference tech- 


niques to predict future reliability in the light of previous 
failure data. 


4. Model Unification. 


By adopting a Bayesian approach, it turns out that one can 
unify models 1, 2 and 8 - the models by Jelinski & Moranda, 
Littlewood & Verall and Goel & Okumoto respectively - under 
a general framework. Observe that this also provides a link 
between the two types of models, since models 1 and 2 are of 
type I whilst model 8 is of type II. 

We begin by recalling the first model, that by Jeliriski & 
Moranda. Each T, is assumed to have a constant failure rate 
A(N-i+1). It is well known that this implies each T, must 
therefore be exponentially distributed with mean (A(N-i+ 1)". 
Now assume that A and/or N is unknown; in true Bayesian 
fashion prior distributions are placed upon them. 

To obtain model 8 by Goel & Okomuto, we let A be 
degenerate at A and N have a Poisson distribution with mean 
8. One can calculate M(t) using the T,’s as defined by Jelinski 
& Moranda, and then by averaging out over N one finds that 
M(t) is indeed a Poisson process with mean: 


y(t) = 0 (1-e™) 


which is the form of u(t) for Goel & Okumoto’s model. 

One can also obtain model 2 by assuming N to be degen- 
erate and to have a gamma distribution. The derivations 
which lead to the above are complex; readers are referred to 
Landberg and Singpurwalla (1985) for the details. 
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5. An Application: Optimal 
Testing of Software. 


The failure models that have been reviewed in the preced- 
ing sections can be used for more than inference or the predic- 
tion of software failure. They can also be applied in the 
framework of decision theory to solve decision problems. An 
important example of such a problem is the optimal time to 
test software before releasing it. This involves the balancing 
of the costs of testing and the risk of software obsolescence 
with the cost of in-service failure, should a bug not be cor- 
rected during the testing period. The following is taken from 
Singpurwalla (1991), in which a strongly Bayesian approach 
is taken. 

To implement a decision theoretic procedure requires two 
key ingredients. The first is a probability model, and here we 
take a generalization of the Jelinski & Moranda model. The 
second is a consideration of the costs and benefits, or utilities, 
associated with a particular decision i.e. the costs of testing, 
the benefits and costs of fixing a bug etc. Decision theory states 
that the optimal decision (in this case time of test) is that which 
maximizes expected utility. 

If the software is to be tested for some time, say T units, 
and then released the problem is to find a T that maximizes 
expected utility. This is called single stage testing. There is a 
more complex, yet realistic, scenario called two stage testing. 
Here the software is tested for T units of time, and then 
depending on how many failures M(T) were observed during 
that test, a decision is made on whether to continue testing for 
a further T* units. The problem here is to find the optimal T 
and T*, with T* to be determined before M(T) is observed. 
Finally there is a third testing scenario, namely sequential 
testing. Here T* is determined after M(T) is observed; this 
procedure can continue for several stages, with T** being 
determined after M(T*) is observed and so on. Here we 
consider the case of single stage testing. Figure 3 is a graph of 
the decision process associated with single stage testing. 

The model chosen in this paper is an extension to Jelinski 
& Moranda’s model. We have 


£. (QINA)=A(N-i+ Ne*N“*Y", 120 


In the previous section we placed prior distributions on 
one of N or A. Now we place priors on both the parameters, 
and say that N has a Poisson distribution with mean 6, A has 
a gamma distribution with scale 1 and shape @ and that N and 
A are independent. 

We now turn to the choice of utility function. The follow- 
ing assumptions are made: 

i) The utility of a program that encounters j bugs during 
its operation is a, + a, ¢*y. 

ii) The cost of fixing a bug is some constant C, 
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Figure 4. 


Time of testing versus expected utility for the software testing 
model. 
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iii) Let f(T) be the cost of testing and lost opportunity to 
time t; here we say f(T) = dT*. 

Note from i) that the utility of a bug-free program is a,+a,, 
and the utility of a program with a very large number of bugs 
is near a,, so that typically a, is a large negative number 
(because there is a great loss associated with software that is 
constantly failing in the marketplace) and a, > 0. Combining 
these assumptions gives us the utility of a program that is 
tested for T units of time, during which M(T) bugs are found 
and corrected, and then released where j bugs are encountered 
by the customer as 


U (TM (1),j) =" x [a, +a, e&*? — C,M (7) - dT") 


where e”" is some devaluating factor. 

Now the two parts of the decision process - the probability 
model and the utility function - are broughi together. We wish 
to find the time f that maximizes expected utility. In other 
words find T such that E(U(T, M(T), j)) is a maximum, where 
we take expectation, using our failure model, with respect to 
M(T) and j. This maximization is quite complex, and must be 
done numerically via computer. The details are found in the 
paper, but the end result is best displayed as a graph of time 
against expected utility (figure 4); in this case one can see that 
the time one should test the software for is about 3.5 units. 





6. Conclusion. 


This paper has attempted to review the main methods, and 
some of the more well known models, that have been used by 
the statistics community in the area of software reliability. The 
first models were almost always based on looking at the failure 
rate of the software; later on the idea of modeling number of 
failures by a Poisson process was used and then most recently 
auto-regressive processes have been suggested as an alterna- 
tive to the failure rate method. Application of the failure 
models, such as to the optimal testing decision problem, is 
another important aspect to the field. 

Earlier it was pointed out that there is almost no data on 
the reliability of commercial software, due to the sensitive 
nature of that information. A possible method of overcoming 
this problem would be to have more interaction between the 
Statistics and computer science communities. In the future, 
such interaction seems essential if models are to become more 
realistic and useful, and it is perhaps surprising that there are 
so few links between the two groups today. 

There still remains much to be researched in this field. In 
the case of optimal testing, plans for two-stage and sequential 
testing need to be developed, whilst the verification of current 
and future models is likely to remain a problem. Nevertheless, 
because of the increasing presence of computers in all aspects 
of our daily lives, the topic of software reliability can only 
become more important in the future. 
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